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Abstract 

Molecular dynamics study on the granular bifurcation in a simple model is 
presented. The model consists of hard disks, which undergo inelastic colli- 
sions; the system is under the uniform external gravity and is driven by the 
heat bath. The competition between the two effects, namely, the gravita- 
tional force and the heat bath, is carefully studied. We found that the system 
shows three phases, namely, the condensed phase, locally fluidized phase, and 
granular turbulent phase, upon increasing the external control parameter. We 
conclude that the transition from the condensed phase to the locally fluidized 
phase is distinguished by the existence of fluidized holes, and the transition 
from the locally fluidized phase to the granular turbulent phase is under- 
stood by the destabilization transition of the fluidized holes due to mutual 
interference. 
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I. INTRODUCTION 



The dynamics of fluidized granular systems has attracted much attention in physics 
communities as a non-equilibrium statistical system [|l]-|3|]. Due to the fact that an elemental 
particle in granular systems is already macroscopic, there are two major differences in their 
dynamics from that of an ordinary molecular system. First, thermal fluctuation does not 
play any role because relevant energy scales for the kinetic and potential energy are much 
larger than the thermal energy. Second, the dynamics is dissipative because the degrees of 
freedom we are dealing with are coupled with microscopic processes that are not treated 
explicitly. Because of these properties, granular systems require an energy source in order 
to be in a steady state and the external gravitational force plays an important role in their 
dynamics. 

In most experimental situations, the systems are excited by a vibrating plate. This is 
called "Granular Vibrated Bed" , in which the effect of gravity is very large @-§f . Granular 
vibrated bed has been studied in many simulation works The convection and 



a linear stability analysis of surface wave [|TT[] in the vibrated bed have been theoretically 
studied. In the conventional setting of vibrated bed, the frequency of external driving 
plays an important role and this gives variety of patterns phases. Experiments on pattern 
formation in granular vibrated bed, have demonstrated the existence of subharmonic stripe, 



square, and hexagonal patterns fl2fl , as well as localized structures called oscillons ||13|| , when 
the frequency and amplitude of the vibration are varied. Comparison between kinetic theory 
and simulation for pattern formation has been done | 14j| . The patterns are similar in spatial 
structure to those observed in fluid dynamical systems, most notably parallel convection 
rolls in a thin layer of fluid heated from below (Rayleigh-Benard convection) and standing 
surface waves in a vertically vibrated liquid layer (Faraday instability). 

However, the vibration complicates the situation. Since energy input from the vibration 
depends on both the amplitude a and the angular frequency u>, the accurate quantities of 
input energy cannot be estimated from the external control parameter T = au 2 /g, which is 
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usually used as the control parameter in the granular vibrated bed. Therefore, it is difficult 
to construct the theory for the dissipative particles system in the non-equilibrium steady 
state. From the numerical point of view, the Distinct Element Method (DEM) ||15[| , which 
is the soft-core particles model, has been often used for the studies of granular vibrated bed. 
However, if we want to obtain the data for the statistical properties during the simulation, 
the sampling data must be taken at the same phase of external vibration cycle. Moreover, in 
the dynamics of dissipative particles, the effect of different time step is not negligible. This 
is a crucial effect to the particle behavior especially when the system is in the high density. 
Therefore, DEM seems to have the practical disadvantages to study the correct statistical 
properties and dynamical behavior in the dissipative system, and is not appropriate as the 
simplest model for the study of the non-equilibrium steady state. The granular vibrated 
bed has these disadvantages to study the statistical characters of the granular system from 
the aspects of both theory and simulation. 

Although there already have been many studies on the granular system for fluidization, 
most of these studies were a collection of the various phenomena from the experimental and 
numerical results, and these were discussed individually for each study. Therefore, it becomes 
difficult to capture what the intrinsic characters of the granular system are. Since the 
study of granular system requires new theoretical ideas beyond ordinary standard statistical 
mechanics or hydrodynamics, a simple ideal model system to compare the theoretical studies 
for the non-equilibrium steady state is required for the progress. Recently, a series of works 
to find some universal frameworks have been appeared (See, Sec. [II A|) . In their study, the 
granular particles are regarded as a collection of inelastic hard spheres that interact with 
each other via hard sphere potential. The assumption of an ideally hard sphere potential 
is also used in kinetic theory and the Boltzmann or Enskog approaches, which facilitates 
comparisons between theory and simulation in the elastic limit. From the numerical point 
of view, even if the typical distance between the particles is short, the molecular dynamics 
event-driven simulation, which is used in the hard-spheres system, is quite efficient and has 
smaller numerical error compared with time-step-driven method (e.g., DEM). 
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Since the gravitational force should have a large effect on clustering, it would be of interest 
to observe how the system behavior changes as the gravity sets in. The main purpose of 
our study is to examine the statistical and dynamical characters of the inelastic hard disk 
system in the non-equilibrium steady state under gravity, and to provide extensive numerical 
results for theoretical studies to construct a base of the universal framework. In this paper, 
we particularly focus on the competition between the effect of external driving and that of 
the gravitational field using the event-driven molecular dynamics simulation, systematically. 

In Sec. [II A| , we will present a short review of a series of previous works using the inelastic 
hard sphere model with the heat bath in the non-equilibrium steady states. Our model is 



described in Sec. [II B| , and we clarify what the advantages and disadvantages of our model 
are in Sec. [II C| and Sec. |II D| , respectively. In Sec. |T|, results of numerical simulations to 
our model are presented. Discussion on the characteristics of the bifurcations are described 
in Sec. |V| Finally, we summarize our work in Sec. [V|. 



II. MODEL 

A. Non-equilibrium steady states 

To achieve the Non-Equilibrium Steady State (NESS) in the granular system, energy 
must be injected into the system from outside. Previous studies for NESS in the inelastic 
hard sphere (IHS) system may be categorized by the following six groups. Firstly, we 
categorize the models by their dimensionality (ID or 2D). Secondly, we consider the way in 
which the system is linked to the heat bath (i.e., heated from the boundaries/bulk). Finally, 
we classify the existence or non-existence of gravity field. We will briefly review these models 
for each category. 
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1. Id model heated from the boundaries without gravity 



Du et al. [16] firstly studied the dynamics of quasi-elastic particles constrained to move 
on a line with energy input from boundaries. They showed that the system violates equi- 
partition of energy and the simple hydrodynamic approach fails to give a correct description 
of the system. Similar models in NESS were also studied from the theoretical point of view 

EM- 

2. Id model heated from the bulk without gravity 

Some papers investigated the properties of the IHS fluid that is heated uniformly by 
random force so that it reaches a spatially homogeneous steady state. This way of forcing 



was proposed by Williams and MacKintosh [21] for inelastic particles moving on a line. Giese 



and Zippelius |22j regarded the internal degree of vibrating rods as a thermal bath. They 
calculated the restitution coefficient r as a function of the relative length of the colliding 
rods, the center of mass velocity, and the degree of excitation of the internal vibrations. 
Puglisi et al. [^3] examined a model consisting of N particles on a ring of length L with 



the Langevin force. They found two regimes: when the typical relaxation time r p of the 
driving Brownian process is small compared with the mean collision time rf the spatial 
density is nearly homogeneous and the probability distribution function (PDF) of velocity is 
Gaussian, but in the opposite limit r p ^> r c one has strong spatial clustering, with a fractal 
distribution of particles, and the velocity PDF strongly deviates from the Gaussian. 

3. Id model heated from the boundaries under gravity 



Kurtze and Hong |24| studied the effect of dissipation on the density profile of ID gran- 
ular gas under gravity. Perturbation analysis of the Boltzmann equation reveals that the 
correction of the density profile along the height y due to dissipation resulting from inelastic 
collisions is positive for < y < y c and negative for y > y c . They also obtained the expres- 
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sion for y c . Using Boltzmann equation and molecular-dynamics simulations, Ramirez and 



Cordero |25| studied the ID column of TV inelastic point particles, in the quasielastic limit, 
under gravity. They showed that the physical properties (e.g., density and temperature 
profile) in the hydrodynamic limit (N — > oo) depend solely on the product (1 — r)N, where 
r is the restitution coefficient, which was confirmed by both theory and simulation. 

4- 2d model heated from the boundaries without gravity 



Esipov and Poschel |26[ studied the kinetic energy distribution function satisfying the 
Boltzmann equation and construct the phase diagram of granular systems on the 2D par- 
ticle system in a circle wall at constant temperature. Their study has suggested that the 
granular hydrodynamics is valid for simple flows or for any flows in dilute and sufficiently 
elastic systems. Grossman et al. [p7|] studied the steady state theoretically. Their predic- 



tions compared well with numerical simulations in the nearly elastic limit. Although the 
system is no longer close to equilibrium even in the low-density limit, the scaling behavior 
of the velocity distributions showed that a hydrodynamic treatment could still be useful. 
The steady state of a low-density IHS gas confined between two parallel walls with same 
temperature is studied by Brey and Cubero f2~8| . Soto et al. |29| also has studied the heat 



conduction between two parallel plates by molecular dynamics simulation of the IHS model. 
Their results show that Fourier's law is not valid and a new term proportional to the density 
gradient must be added. 

5. 2d model heated from the bulk without gravity 



Peng and Ohta |30| studied the 2D granular system where an external driving force is 
applied to each particle in the system in such a way that the system is driven into a steady 
state by balancing the energy input and the dissipation due to inelastic collisions between 
particles. Experimentally, steady state of clustering, long-range order, and inelastic collapse 
are observed in vertically shaken granular monolayer, which is correspond to the system 



heated from the bulk in the limit of the frequency ui — > oo pi|. In 2D the model may be 



considered to describe the dynamics of disks moving on an air table. A similar IHS model 



with random external accelerations has been studied by Bizon et al. [32] to test continuum 



theories for vertically vibrated layers of granular material. To investigate the long-scale 



correlation and structure factor, van Noije et al. [BSJ present a theory of the randomly driven 



IHS fluid in 2D simulation, and characterized its non-equilibrium steady state. As another 



way to realize the steady state, Komatsu |34[ introduced a new observation frame with 



a rescaled time. This operation is the same as the velocity scaling method, which is often 
used as the temperature control method of the molecular dynamics simulation in equilibrium. 
He found that a phase transition with symmetry breaking, which is the circulation of the 
particles in the box with elastic walls, occurs when the magnitude of dissipation is greater 
than a critical value. 

6. 2d model heated from the boundaries under gravity 



In 2D system under the gravity, Luding et al. |35j compared the results of DEM and IHS 
model. The system of their study is driven by vibrations like a granular vibrated bed. The 
results of DEM simulation were good agreement with that of event-driven simulation in the 
fluidized regime. They also found that the results of DEM simulation depend strongly on 
the time t c during which the particles are in contact. Although there exist only a few recent 
laboratory experiments for 2D system under the gravity, an experiment on the stainless 



steel particle system in two dimensions was studied by Kudrolli et al. [p6 |. The system was 
driven by the periodically moving wall and slightly inclined to apply the weak gravity. This 
experiment showed that the clustering occurs around the opposite side of driving wall when 
there is no gravity, but the cluster migrates downward when the system is inclined even by 
a very small angle. 
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B. IHSHG model 



To study the statistical properties of the granular system under an external field with 
energy input in NESS, we propose a simple model system on the 2D inelastic hard disk 



system heated from the boundaries under the gravity [37]. 

The system considered here consists of inelastic hard disks of mass m and diameter d 
in 2D space under uniform gravity. We employ the periodic boundary condition in the 
horizontal direction. The particles collide with each other with a restitution coefficient r. 
All disks are identical, namely the system is monodispersed. For simplicity, we neglect the 
rotational degree of freedom. Because of the hard-core interaction, collision is instantaneous 
and only binary collisions occur. When two disks, i and j, with respective velocities and 
Vj collide, the velocities after the collision, and vj, are given by 

v^ = v < -^(l + r)[n.(v i -v i )]n (2.1) 
v' j = v j + ±(l + r)[n.(v i -v j )]n, (2.2) 

where n is the unit vector parallel to the relative position of the two colliding particles 
in contact. Between the colliding events, the particles undergo free fall motion with the 
gravitational acceleration g = (0, —g) following parabolic trajectories. The system is driven 
by a heat bath with temperature T w at the bottom of the system; a disk hitting the bottom 
bounces back with velocity v = (v x , v y ) chosen randomly from the probability distributions 
(p x (v x ) and <j>y{v y ) @; 



j ™ TtViP' 

Mvx) = J , ex P(~ , Z ) (-°° <v x <oo) (2.3) 
V 27rk B T w 2k B T w 
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(j) y (v y ) = — -v y exp(- y ) (0<u y <oo), (2.4) 

KbJ-w ZKbJ-w 

where ks is the Boltzmann constant. Note that the wall temperature T w here is just a 
parameter to characterize the external driving and is not related to the thermal fluctuation. 
We call this system "Inelastic Hard disk System with a Heat bath under weak Gravity" 
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(IHSHG) [[37]] to emphasize the importance of the competition between excitation by the 
heat bath and the gravity. 



C. Advantages of IHSHG model 

The IHSHG model is so simple that the system is completely characterized with only four 
dimensionless parameters: the restitution coefficient r, the driving intensity A = ksTw/mgd, 
the system width in the unit of disk diameter N w = L x /d, and the particle number of the 
initial layer thickness = N/N w , where L x is the system width and N is the total number 
of disks. 

Now we compare our IHSHG model with other typical dissipative system in NESS. 
Rayleigh-Benard(RB) convection is the typical continuous fluid system, where the energy 
source is the heat bath at the bottom, the dynamics is described by Navier-Stokes equation, 
and the system is driven by the difference of temperature between the top and the bottom. 
The dissipation comes from the viscosity of fluid. In this system, the control parameters are 
well known as Rayleigh Number and Prandtl Number. There are some previous MD studies 
of RB convection using elastic hard disk |39| , ^0| . 

The ordinary granular vibrated beds by using DEM |7|-f| are the discrete soft particles 
systems; the energy source is the vibration of the bottom, and the dynamics obeys Newton's 
equation of motion under the gravitational field. The dissipation appears between inelastic 
particles, and the control parameter is Y = au 2 /g, where a and lo are the amplitude and the 
angular frequency respectively. Recently, Aoki and Akiyama [^IJ proposed that the state of 



granular convection involving multipairs of convection rolls is governed by a dimensionless 
parameter A = = ^=p, where D and a are the self-diffusion constant and the energy 
diffusivity a, respectively. 

Our IHSHG model is the discrete hard disks system, the energy source is the heat bath 
at the bottom and the dynamics obeys Newton's equation of motion under the gravitational 
field. The dissipation comes from inelastic collision of particles, and the control parameter 
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is A = ksTw/mgd. The present system is analogous to the ordinary granular vibrated bed 
but simpler because it does not have an external time scale. Since IHSHG model has no 
external periodic cycles and there are only four dimensionless parameters in the system, 
IHSHG model is easy to deal with numerical analysis and theoretical approach. 

As mentioned briefly in Sec. [11 A H, there is a recent experiment [BH] on the stainless 



steel particle system in 2D similar to IHSHG model. Based on these experimental results, 
we firstly reproduced the clustering phenomenon and the behaviors of NESS to simulate 
IHSHG model. Consequently, we found that IHSHG model could reproduce Kudorolli's 
experiment qualitatively. Note that IHSHG model might be categorized in 2D model heated 



from the boundaries under gravity described in Sec. [II A 6 



D. Some problems of IHSHG model 

Although it is simple and easy to deal with, there are some problems in IHSHG model. 
Since IHSHG model is based on the inelastic hard-sphere model, all problems come from a 
reason that the potential between two colliding particles is unusually stiff. The instantaneous 
collisions imply that an interaction is zero when the particles are not in contact and suddenly 
becomes infinite when the distance between the colliding particles is zero. Therefore, the 
momentum exchange occurs instantaneously. However, in a real system, the situation is 
different, each contact takes a finite length of time and the force is large but finite. The 
infinitely stiff hard-sphere interaction is an idealization or simplification of smooth repulsive 
pair-potential. 

When the restitution coefficient is low (i.e., a dissipation rate is high) or kinetic energy 
is relatively smaller than potential energy under external gravity force, the number of col- 
lisions during a finite time often diverges. This is called "Inelastic Collapse" . It was firstly 
discovered while studying the ID model system of dissipative particles with a elastic wall 
?2|,|43[], thereafter in 2D |44]]. Although the inelastic collapse is thought as artificial pathol- 



ogy, several proposals exist to avoid it |22|H3-BU[ . The velocity dependence of the restitution 
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coefficient was experimentally measured for binary collisions, and the quantitative study of 



the internal modes of a granular particle is under way [SI]. However, the inelastic collapse 
has never appeared when the dissipation is low and the system is highly excited. Therefore, 
IHSHG model for the parameter region studied here does not need to take care of it. 

Since the dynamics proceeds by only the event of binary collision, multiparticle inter- 
actions are impossible to realize. In a real system, each contact takes a finite time so that 
multiparticle contacts are possible. Nevertheless, in the fluidization study, multiparticle 
collision can be thought as a rare event. Hence, we can ignore this effect. 

Another problem with the inelastic hard-sphere model is that the static limit does not 
exist. For example, in the framework of the inelastic hard-sphere model, a particle cannot 
rest motionless on the ground. In the realistic granular systems, the potential and elastic 
energies will approach constant values while the kinetic energy tends to zero. If the particles 
would be dissipative hard particles, the inelastic collapse will occur before the kinetic energy 
vanishes. The contact forces and stresses are not properly defined so that the system will 
never reach a static configuration. However, since IHSHG model has a heat bath at the 
bottom of the system, the kinetic energy always stays finite. In the simulation of NESS, the 
dynamics has little effect from this problem. 



III. NUMERICAL SIMULATIONS 



In this section, we present the results of numerical simulation on IHSHG model using 



the simple and efficient event-driven algorithm |52| , which allowed us to study the statistical 



and dynamical properties for a wide range of parameters. 



A. Numerical setting 

Most of the simulations were performed on the system with r = 0.9, N w = 200, and 
Nh = 20 (i.e., N = 4,000). The driving intensity A = ksT/mgd was varied from 25.0 to 
909.1 to study the change of the system behavior. In the actual simulation, we reduced 
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any variables to dimensionless ones by the following three basic units (ksT, m, d), which 
are all fixed at unity. Note that unit time r(r = yf^fO is also unity in these simulations. 
To ensure that the system had reached a steady state, we relaxed the system until the 
total energy did not drift. Therefore, the system reaches the non-equilibrium steady state. 
Thereafter, we started the simulations with the length of 100,000 collisions per particle; we 
could obtain the data with small statistical error. In addition, we performed the simulations 
with N w = 400 (i.e., N = 8,000) and N w = 1,000 (i.e., N = 20,000) to study the size- 
dependence of the system behavior. These run lengths were 30,000 and 10,000 collisions per 
particle, respectively. The discussion for the size-dependence of the system is presented in 
Sec. 0. 

B. System behavior upon changing external driving 

1. Packing fraction and excitation ratio 

In order to characterize the steady states of the system, firstly, we measure the maximum 
packing fraction A max at the height y = HAmax as a function of the driving intensity A. In 
the system with (r,N w ,N h ) = (0.9,200,20), there are two cusps around A ~ 200 and 
380(o in Fig. [I]), suggesting phase transitions with changing A. These transitions should be 
related to the excitation structure of the state, therefore, we define the excitation ratio fi 
by n = fC/U, namely, the ratio of the total kinetic energy /C = (m/2) J2i v f to the potential 
energy U = rngJ^iVi- The cusps appear at the same points for fi, which confirms that some 
transitions occur at these points (x in Fig. |l|). 

2. Snapshots 

In order to conceptualize the underlying physical mechanism of the system behavior in 
each phase separated by the transitions, three snapshots for typical values of A for each 
phase are shown in Fig. |2], in which interval time between them are all lOOr. 
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For A = 181.8 (al — a3 in Fig. |2|), most of the particles are aggregated around the bottom 
of the system with an almost close packed density, and the state of a dense layer is relatively 
stable, which means the potential energy is dominant and the system is in weakly excited 
states. We call the phase in A < 188.7, the condensed phase (CP). It can be seen from the 
snapshots, however, that collective motion appears on the surface of the layer. 

In the second phase, A = 250.0 (bl — b3 in Fig. |2|), the dense layer packing is locally 
broken by excitation of the heat bath. The high-speed particles are blown upward from the 
holes in the layer. In Fig. |2], we see only one hole in the system, but for a larger system, 
there are certain cases where we can observe more than one hole. The holes migrate and 
occasionally become more active; sometimes they almost close temporarily, but the structure 
is fairly stable. We call this phase the locally ffuidized phase (LFP). 

For the case of A = 666.7 in the third phase (cl — c3 in Fig. |2|), the layer is completely 
destroyed. The average density is quite low, but it is very different from the ordinary 
molecular gas phase. The density fluctuation is large and this fluctuation causes turbulent 
motion driven by the gravity. At some time, the whole system gets excited with some 
relatively smaller density fluctuations, but the very next moment, a large proportion of the 
particles travels downward and forms a layer like structure. This structure, however, is 
destroyed immediately. We call this phase the granular turbulent phase (GTP). 

3. Collision rate 

Since the simulation was performed in the non-equilibrium steady state, the collision rate 
takes constant value during the whole run. Therefore, when the total number of collision 
C total {= 100,000 x N) was fixed, the total time T to tai for each run were different for each 
A. We found that the system in the non-equilibrium steady state can be also characterized 
by the collision rate because it seems to correspond to the pressure of molecular system 
in the equilibrium. Though in the elastic hard disk system the collision rate is only the 
function of packing fraction of the system, in the inelastic hard disk system it depend on 
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the competition between the inelasticity r and the external driving A. 

In Fig. [| the collision rate (C tata i/T tota i) vs. the external driving A in the log-linear 
scale is shown. We found that the collision rate decrease exponentially in LFP and GTP 
(A > 188.7) with different exponents. However, the collision rate deviates from exponential 
fitting in CP. It diverges when the system is in the weakly excited state (A — > 0). 

4- Probability distribution function of velocity 

These inhomogeneous behaviors should result in a non-Maxwell-Boltzmann distribution 
of velocity. Figure |] shows the horizontal velocity Probability Distribution Functions (PDF) 
at the horizontal layer around the height of the maximum packing fraction (y = HA max ) 
for each phase in the log-linear scale. The data are plotted in the unit of ^k^T /m. The 
distributions for A = 250.0 and 666.7 deviate from the Gaussian and are more or less in 
exponential form in the tail region. It can be seen that the distribution for A = 250.0 in 
LFP consists of two parts; the central part that originates from the dense layer packing, and 
the wide tail that comes from the fluidized holes. The central parts for all cases are very 
close to Gaussian (inset of Fig. |j). 

5. Flatness parameter 

In order to quantify the deviation from Gaussian distribution, we next calculated the 
flatness parameter / as a function of the height y defined by 

M = «)/{vl)\ (3.1) 

which is often used in the study of the intermittency of the fluid turbulence. The dependence 
of the external driving of the maximum flatness f{Hf max ) at the horizontal layers with the 
thickness of disk diameter d is shown in Fig. |5]. Over the whole region of external driving 
A, the value of f{Hf max ) is different from 3, which correspond to the value for the Gaussian 
distribution, but it is remarkable that / becomes very large, as large as 20, in LFP. In Fig. ||, 
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the statistical error bars are also plotted for each data. Note that the statistical errors of 
the time averages for molecular dynamics simulation were estimated by the block average 
method described in ref. ||53|| . 

6. Fitting functions of velocity PDF 

The unexpected large value of /, however, can be understood as follows |3"7fl . Assume 
the velocity PDF 0(t>) has two components, the narrow Gaussian distribution 4>g( v ) with 
the weight 1—p and the broader stretched exponential distribution <f>s(v) with the weight p; 

<j>(v) = p<t> s {v) + (1 - p)<P G (v), (0 < p < 1), (3.2) 

where <fis( x ) an d 4>g{ x ) are defined by <ps ( x ) = (/3/(2xoT(l//3)))exp(—\x/xo[) and 4>g{ x ) = 
(1/(v / 27tctg)) exp (— x 2 /(2a G )), respectively. Then, the flatness of this distribution is given 
by 

[T(5/(3)/T(l/f3)}xlp + 3(l-p) „^xo 
{[T{?>/(3)/T{l/(3)]xlp + {l-p)f' ° °g 

because the second and the fourth moment of (j>s{x) are given by r(3//3)/T(l//3) • Xq and 

r(5//3)/r(l//3)-Xo, respectively. If we fit the parameters in Eq. ( |3.2p to the data for A = 250.0 

for example, we obtain oq = 0.048, xo = 1.7, (3 = 0.85, and p = 0.25, which yields / ~ 19.0 

(Fig. U). The enhancement of the flatness by the superposition of the broader distribution 

can even be drastic for a small value of p as can be seen if Eq. ( |3.3|) is plotted as a function of 

p. This observation indicates that / can be used as a sensitive index to detect the appearance 

of a small weight of the broad component in the distribution. The sharp rise of / in Fig. |5| 

around A ~ 200 is clear evidence of the appearance of the fluidized holes. 

The other PDFs for each A in LFP and GTP are also fitted by changing p (e.g., A = 476.2 

in Fig. |6|). From these p, we obtain the plot of p for each A in Fig. [?|. Near CP-LFP critical 

point Aql, the weight of p takes a certain finite value, not (Fig. 0). In LFP, the weight p 

increases linearly when A increases, which means the fluidized holes region increases linearly. 

However, there is a plateau p ~ 0.5 when the transition from LFP to GTP occurs. 
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C. Height dependencies of flatness and packing fraction 



The height dependencies of the packing fraction A(y) and the flatness f(y) are calculated 
in the horizontal layers with the thickness of the disk diameter d. In Fig. ^, we plot the A 
dependencies of Hf max and HAmax, which are the heights where the flatness and the packing 
fraction take their maximum values. In LFP and GTP, the both of the parameters reach their 
maximum values at the almost same height (Hf max ~ HAmax), but in CP the flatness shows 
its maximum value at the lower position than the packing fraction does (Hf max < HAmax)- 
In CP, the value of the packing fraction where the flatness takes its maximum value is very 
close to the value of the Alder transition point A c (~ 0.7) in the elastic hard disk system 
(Fig. P). The fact is somewhat intriguing, but we could not find obvious reason for it. 

In the simulation of the granular vibrated bed, Taguchi and Takayasu |54| showed that 
velocity PDF becomes power distribution in the fluidized region, but Gaussian in the solid 
(close packed) region. On the other hand, Murayama and Sano |55j found that the flatness 
varied from Gaussian (/ = 3) to exponential (/ = 6) when the packing fraction increases. 
These results seem to disagree with each other. However, in CP of IHSHG model, we 
obtained two different characteristics about velocity PDF; (i) The region < y < Hf max , 
where the packing fraction increases from to A c and the flatness also increases, (ii) the 
region Hf max < y < HAmax, where the packing fraction (> A c ) reaches nearly the closed 
packing and the flatness decreases. 



D. Spatial velocity correlation function 

In the region (ii) of CP (i.e., Hf max < y < HAmax), we found the Spatial Velocity 
Correlation Function (SVCF) C VxVx (R) = (v x (x + R)v x (x)) takes positive value for all relative 
displacements R(0 < R < L x /2) (e.g., A = 175.4 (CP) in Fig. HQ). Thus, the particles flow 
collectively in the same horizontal direction. 

In Fig. [rg, SVCFs at the height H Am ax in LFP and GTP are also shown. In LFP, the 
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negative minimum of SVCF becomes larger when A is increased, which means that the size 
of fruidized hole grows large. In GTP (i.e., A > 357.1), there is no negative minimum in 
C VxVx (R) for all R less than half of the system size(L x /2). This means the size of fruidized 
hole grows larger than the half of system size. 



E. Displacement of particles 



In Fig. [XX], the time-dependence of displacement vectors for each particle from the initial 
position are shown at t/r = 100,200,300,400 when A = 25.0. The displacement vectors 
from the dark region indicate that the direction of the displacement is right. This is a 
clear evidence of the existence of the domain of collective particle motion. The feedback of 
collective motion can clearly be seen during the time t/r = 200 ~ t/r — 400. 

In Fig. |T2|(a), we also show a snapshot of the displacement vectors for each particle at 
t/r = 150 in A = 175.4. The domain size of collective motion becomes smaller than that 
of Fig. [11]. Since the particle has a room to move freely compared with a weakly exited 
state (e.g., A = 25.0), the convective motion appears near the critical point Acl- This 
convectional mode is similar to the granular vibrated bed. 

In LFP (Fig. |TJ(b)), particles at the region of broken layers behave such as vortex-like 
stable convection. The most of excited particles drops on the condensed layer near the 
broken region. 

In GTP (Fig. pT2](c) ) , the particles at the broken layer are highly excited, which can reach 
an another side of the condensed layer through the periodic boundary condition. This motion 
must influence the results on the velocity PDF and spatial velocity correlation function. 



F. Horizontal averaged particle current in CP 

Since the excluded volume effect is dominant in CP, the collective motion of particle 



appears. Figure 11 shows that the displacement of particles in CP move in the same direction 
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during the short time scale. We can also confirm the above result quantitatively that the 
averaged particle current field summed up in the packed region (ii), 

jfo,Af) = -i- £ 5>(At)5(x,(At)-x 9 ), (3.4) 

s region(ii) i 

where N s is the sample number of particles, x g is grid with the width of the disk diameter 
d and At = N j ' [Colli sionRate). Figure [T3] shows that the time-evolution of the horizontal 
averaged particle current for one particle in A = 25.0 and A = 175.4. In a relatively short 
time scale, the particle currents in the packed region take positive value. However, we found 
the currents fluctuate around zero in a long time scale. The feedback time scale of the 
current becomes longer when A is near the critical point. Therefore, the collective motion 
appears strongly when the system is in the weakly excited state. 

G. Surface wave in CP 

In CP, the dense layer packing is formed and it is stable in time because excitation is 
weak and the potential energy is dominant. However, it is apparent that collective motion 
appears on the surface of the condensed layer. What is the mechanism of the surface wave- 
like motion in CP? There must be no surface tension in the macroscopic dissipative system, 
therefore the restoring force should originate from the balance between the gravity and the 
excitation, but it is not clear if it could be defined as an ordinary gravitational wave in fluid. 

To study the dynamical behavior of surface wave in detail, simulations were done on 
the driving intensity A varied from 25.0 to 188.6. We performed simulations until the time 
T = 10240r for each A. Figures |T4"|(a) and (b) show the spatio-temporal pattern of surface 
height of the condensed layer in CP far from the CP-LFP critical point ((a) A = 25.0 and (b) 
100.0). And Fig. [TJ(c) shows the pattern near the CP-LFP critical point ((c) A = 188.6). 
The horizontal and vertical axes correspond to the space expansion (x = ~ L x ) and the 
time evolution (T = ~ 2000r), respectively. There are many domains in the spatio- 
temporal patterns, which represent the collective motion of surface appears. The typical 
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domain sizes, which are estimated by Fig. |14Kc), are about Ax ~ L x /2 and At ~ 850r. 
These are correspond to {k(= 2ir/(L x /2)),u(= 2tt/850)) = (0.0628,0.00739), respectively. 
To study this collective motion quantitatively, we calculate the dynamical structure factor 
of the surface wave for each run. The dynamical structure factor S(k,u) is described as 



where h y (x,t) is the surface height of condensed layer; N gx (— 200) is the total number 
of grids with the width of disk diameter d in the horizontal direction and the discrete 
positions are taken by Xi = (L x /N gx )/2 + (i — 1)(L X /N gx ), (i = l...N gx ). When A = 188.6, 
the strongest peak of dynamical structure factor in the k — uo space appeared at (k, uS) ~ 
(27r/(L x /2), 0.0065). Thus, we confirmed these values were consistent with the estimation 
of the domain sizes in the spatio-temporal pattern. When A becomes low, the domain size 
in the vertical direction becomes small (Fig. 0). On the other hand, the domain size in the 
horizontal direction is not changed very much. This means that the periodical cycle of surface 
wave becomes shorter in low A. These results are also confirmed by the computation of the 
dynamical structure factor for each A. At a result, when the strongest peak (k max , uo max ) was 
plotted in k — uj space for each A, k max stays around k = 2n/(L x /2) and increases slightly 
upon decreasing A, but u max becomes large as A go to 0. 

We conclude the periodic behavior of the surface wave in CP shows the critical slowing 
down when A — > Aql- Note that the value of u at the critical point takes a certain finite 
value, not 0. It can be understood that the transition occurs before the break down of the 
condensed layers, which means the transition seems to be subcritical. We also found that 
the amplitude of the surface fluctuation becomes large near the transition point Aql- This 
is an interesting point that the system shows large fluctuation around the transition like the 
critical phenomena in equilibrium. 




x exp (—ik(xi — Xj)) exp (iut)dt, 



(3.5) 
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H. Dynamics of fluidized holes 



The localized excitation in LFP should resemble a circle in 3D system and reminds us 



of an oscillon [13|, but they are different; the external vibration frequency is essential for 
the oscillon dynamics, but the localized excitation here does not have such a characteristic 
frequency. In the large N w simulation, we can observe the merging process of two or more 
excitations, but their mode of interaction is not clear. An interesting question here is if 
the transition from LFP to CP and/or the transition to GTP can be understood in terms 
of the local excitation; if LFP transform to CP when the distance between the excitations 
diverges? Does LFP transform to GTP at the point where they condense? 

To investigate the excitation dynamics, we calculate the density fluctuation for each 
position with the width of particle diameter d in the system. The driving intensity A was 
varied from 200.0 to 476.2, in which the systems are in LFP and GTP. We performed 
simulations until the time T = 10240r for each A. The particle number in each column, 
which is summed up from the bottom to the top, defines the density. If the particle number 
is less than the averaged layer number Nh, we regard the column as the fluidized holes (black 
region in Fig. |15|). If the column has the particles more than the averaged layer number Nh, 
these columns are regarded as condensed layers. This criterion seems not to be a rigorous 
definition of fluidized holes and condensed layers, however, we confirmed that the fluidized 
holes were distinguished very well. 

In Fig. [15|, the spatio-temporal patterns of the density fluctuation are shown in LFP 
(N w = 200; (a) A = 200.0 and (b) 250.0) and GTP (N w = 200; (c) 476.2). Near the 
critical point Aql, the lifetime of a fluidized hole does not continue through the simulation 



(Fig. |15|(a)). When A is increased, the region of a fluidized hole increases (Fig. |15|(b)). 
These results are consistent with the increase of the weight of the broader distribution p 
for the fitting functions of velocity PDF (Sec. p 1 1 B 6|) . Above LFP-GTP transition point 
(Fig. |T5|(c) ) , the local excitation domain fluctuates largely and sometimes expands to the 
length of the system. We found that the spatio-temporal pattern of a fluidized hole changes 
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significantly through the transition from LFP to GTP. 



IV. DISCUSSION 



A. Size dependence 



In this subsection, we discuss the changes of the statistical properties in IHSHG model 
upon increasing the system size N w . 



We performed the almost same simulations described by Sec. |1| for the systems with 
larger N w , which means the aspect ratio in the system is changed. To examine the size 
dependence, the numerical simulations for various A with N w = 400 (i.e., iV = 8,000) and 
N w = 1,000 (i.e., iV = 20,000) are performed. The simulation lengths for these data are 
30,000 (N w = 400) and 10,000 (N w = 1,000) collisions per particle, respectively. We found 
that the N w dependencies of the maximum packing fraction and the excitation ratio do not 
change when N w is increased. However, the N w dependence of the flatness at the height 
HAmax clearly becomes larger in LFP and GTP. 



In Fig. 16, the velocity PDFs <fi(v) with fitting function in N w = 200 and A^ = 400 
are shown when A = 250.0 (LFP). The difference of velocity PDF between N w = 200 and 
N w = 400 is very small. However, the flatness becomes clearly larger when N w is increased 
(e.g., N w = 200, / ~ 19.0; N w = 400, / ~ 23.0). We change (3 for N w = 400(/3 = 0.82) to fit 
slightly broader distribution than that of N w = 200(/3 = 0.85). The other fitting parameters 
for N w = 400 remain the same as in N w = 200, i.e. oq = 0.048, xq = 1.7, and p = 0.25 



(Sec. f l 1 B 6| ) . The fitting function seems good agreement with the simulation data. Using 



these parameters, we calculated the flatness value by using Eq. (|3.3|) . In the inset of Fig. [16|, 
the flatness values both j3 = 0.85 and (3 = 0.82 are plotted as functions of the weight p. 
We found that the flatness values with f3 = 0.82 are larger than that with (3 = 0.85. This 
is consistent with the fact that the flatness in N w = 400 takes larger value than that in 
AL = 200. 
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Although there is one hole in the system when N w = 200 (Sec. |I7l|) , there should be 
more holes in the limit N w — > oo to reach finite density of holes. What is the average 
density of the holes? It might have some relationship to the width of the surface wave in CP 
(Sec. [Ill G| ). In Fig. [I7|, the spatio-temporal pattern of the density fluctuations are shown in 
LFP (N w = 400; (a) A = 200.0, (b) 250.0, and (c) 333.3) and GTP ((d) A = 454.5, (e) 588.2, 
and (f) 833.3), respectively. Near the critical point Acl(A = 200.0), the lifetime of a fluidized 
hole is shorter than the simulation length (Fig. 0(a)), which is similar situation in the case 
of N w = 200. When A is increased, the width of the holes becomes wider. However, the 
system seems to contain two fluidized holes (Fig. |17|(b),(c)), which is the different situation 
if it is compared with the case of N w = 200. In N w = 1000, we confirmed the fact that there 
are multiple holes in the system. As mentioned above, since the maximum packing fraction 
for each A is not dependent on the system size, the ratios of the total hole region to the 
system size for each N w are almost equal. However, the average size of hole is not dependent 
on system size, since it is determined by the driving intensity A. At the result, the systems 
with large N w have several holes. Above the LFP-GTP transition point (A > 357.1), the 
local excitation domain largely fluctuates (Fig. |T7|(d),(e)). We can see the high-speed motion 
of the high density cluster when A = 833.3 (Fig. [Tj](f)). When N w is increased, since there 
are many fluidized holes in the system, the fluctuation of hole width might become larger. 
In this situation, velocity PDFs for each hole have different variances for the width of holes, 
which means various broader distributions of horizontal velocity coexist with in the system. 
In general, the flatness takes large value when the spatial heterogeneity become large. This 
spatial heterogeneity of the velocity PDFs at fluidized holes might be the reason why the 
flatness takes larger value when N w is increased. 

B. CP-LFP transition 

In this subsection, we summarize and discuss the results of numerical simulation about 
the transition from CP to LFP. 
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We found that the value of the maximum flatness at y = Hf max shows a sharp rise around 
the transition point between CP and LFP {Acl ~ 188.7). The flatness / becomes larger 
than 20, and we demonstrated that this unexpected large value of / comes from the fact 
that the velocity PDF in LFP consists of two components; the narrow Gaussian distribution 
that comes from the dense layer and the broad stretched exponential distribution in the 
fluidized holes (Sec. [II [B 6|) . Actually, velocity PDF <fr(v) in LFP can be expressed very well 
by the sum of the narrow Gaussian <pc{v) and the broader stretched exponential 4>s{v) as 
4>{v) = p(f>s(v) + (1 —p)4>g(v), with the weight p (0 < p < 1). Therefore, CP and LFP are 
distinguished by the existence of the fluidized holes in the latter. 

Near the critical point Acl i n LFP (188.7 < A < 204.1), we found that the lifetime of a 
fluidized hole is shorter than the simulation length (Fig. |l"5|(a) and Fig. |17](a)). On the other 
hand, for 204.1 < A < 357.1, fluidized holes persist in time. In LFP, the majority of the 
particles are in the condensed layer, and it is punctured by the fluidized holes, from which the 
high-speed particles are blown. We observed the stable flow of particles between the fluidized 
holes and the edges of the condensed layer (Fig. |T2|(b)), and this stable flow stabilizes the 
isolated fluidized hole in the condensed layer. The horizontal velocity correlation between 
the opposite edges of the holes becomes negative, and this should be the origin of the 
negative minimum in the spatial velocity correlation function (Fig. [ToD . We observed the 
size of the holes gets larger when A is increased. This result can be understood by the fact 
that the size is determined by the distance that the blown up particles reach, which is, in 
turn, determined by the driving intensity A. 

The time scale of the surface fluctuation in CP becomes larger as A approaches the 
critical point Acl from below. On the other hand, the length scale seems to stays fairly 
constant. From the detailed analysis of the dynamical structure factor of surface wave, the 
typical frequency u does not vanish at the critical point A C l, which suggests the transition 
between CP and LFP is subcritical. 
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C. LFP-GTP transition 



As the driving intensity is increased, the density and the size of the fluidized hole increase, 
and at A^g ~ 357.1 neighboring holes start interfering with each other. At the transition 
point Alg, the blown up particles at one hole barely reach the neighboring holes. As we 
have discussed, the distance that the flight of blown up particles also determines the size of 
the holes, therefore, at A = Alg, we expect the half of the condensed layer structure should 
be destroyed by the fluidized holes. In Fig. [I], we found that the averaged packing fraction 
at A = 357.1 is A ~ 0.48, which is a little larger than the half of that of the value of the 
closed packing fraction (Acp ~ 0.901) in 2D. 

In the system with N w = 200, for which most of the simulations were done, there appears 
only one fluidized hole. In this situation, the interference between holes takes place only 
through the periodic boundary. Therefore, the LFP-GTP transition occurs when both the 
distance and the size of holes become a half of the system size L x /2. And the spatial velocity 
correlation function C VxVx (R) between the opposite edges of the hole becomes negative. In 



Fig [10], we actually confirm the results on the negative value of SVCF C VxVx (R = L x /2) 
when A > A LG . 

Another peculiar aspect of the system with N w = 200 is that the weight of the stretched 
exponential distribution in velocity PDF (i.e., p) stays ~ 0.5 for a finite range of A in GTP 
(Fig. 0). This plateau behavior may be also interpreted as the feedback effect of the periodic 
boundary condition which enforces the number of holes to be one for a certain parameter 
region of A. In the larger systems with N w = 400 and 1,000, we observed more than one 
holes in the system, and found two different holes interfere with each other in GTP. 

The interference between holes destroys the stable particle flow within an isolated flu- 
idized hole (Fig. [L2](c)), and the holes are destabilized dynamically (Fig. [15] and Fig. |T7D . 
In Fig. |15|(c) and Fig. |17](d),(e),(f), which are the spatio-temporal patterns in GTP, we 
found many high-speed moving small clusters. These clusters divide the region of the flu- 
idized holes (i.e., black region) into several parts. In the patterns of LFP (Fig. |l|(b) and 
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Fig. |17|(b),(c)), there are no such moving clusters. Therefore, these high-speed moving small 
clusters are thought as the peculiar phenomena in GTP. When these clusters move across the 
fluidized holes in the spatio-temporal patterns, the lifetime of the fluidized holes becomes 
shorter. Although there is still room for the definition of quantitative values to distin- 
guish the phases between LFP and GTP, we found that the maximum lifetime of fluidized 
holes in the spatio-temporal patterns as a function of A changes drastically between LFP 
and GTP. Thus, we conclude that the LFP-GTP transition is interpreted as the dynamical 
destabilization transition of the fluidized holes. 

The particle dynamics in GTP are very different from that of LFP. The average density 
is quite low, but it is very different from the ordinary molecular gas phase. The density 
fluctuation is very large and this fluctuation causes turbulent motion due to the gravity. 
Even for very large A, the high-density clusters which look like the condensed layer are 
formed temporally. It might be also different from ordinary fluid turbulence. In the finite 
system with finite height, turbulence may disappear when the driving A is large enough, but 
in the system with infinite height, it appears that turbulent motion persists however large 
the driving A is. 

D. Various limits of IHSHG model 

The dynamics of IHSHG model is very different from those of not only an ordinary 
molecular system but also of a conventional granular vibrated bed. Let us discuss the 
relationship between the present model and the vibrated bed, where the system is driven 
by a vibrating plate at a given frequency. In such a case, the amplitude a and the angular 
frequency u of the vibration are two independent parameters, but the control parameter is 
usually taken as the ratio of the accelerations T = au 2 / g. In most experimental situations, 
the external frequency determine the time scale of the system and the collision interval is 
directly related to 1/uj. In the present case, the ratio of the energies A = kBT w /mgd is the 
only control parameter and no external time scale is imposed. This situation may correspond 
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to that in the vibrated bed when 1/u is much shorter than any other of the relevant time 
scales in the system. 

Since we fixed the parameters of IHSHG model at Nh = 20 and r = 0.9, the dependence 
of the behavior by changing Nh or r is not obvious. When the restitution coefficient is fixed at 
unity, we can obtain the density profile of point particles (i.e., the diameter d — > 0) within 
the framework of the statistical mechanics. Under the situation where the temperature 
and the gravity are uniform in the vertical direction, the density profile as a function of 
the height y is exponentially decreasing function ~ exp (—mgy/kBT) and the excitation 
ratio fi is always unity. Since the point particles can be compressed infinitely, the profile 
simply shifts. Therefore, no phase changes occur in the system when the driving intensity is 
changed. However, in the case of finite diameter of particle d, which cannot be compressed 
infinitely, the density profile becomes different from that of point particles when the driving 
intensity becomes relatively low. Under a certain intensity of driving force, a fraction of hard 
particles condenses from the bottom toward the surface. Within the framework of Enskog 
theory [pEfl , the density profile for hard spheres remains finite at the close packed density 
near the bottom and the exciation ratio /i becomes less than unity. It is difficult to predict 
the system behavior when Nh is changed, because Nh dependence of the system might be 
significantly changed for various restitution coefficient r. However, in the case of Nh = 1, 
the system is almost the same situation as ID model heated from the bulk (Sec. [H A 2| ). 
Systematic studies of the system behaviors by changing these parameters will be the subject 
of future work. 

V. CONCLUDING REMARKS 

In this paper, we performed numerical simulations for a two-dimensional inelastic hard 
disk system heated from the boundaries under gravity (IHSHG). Upon increasing the heat 
bath temperature, the system exhibits three distinct phases, namely, the condensed phase, 
the locally fluidized phase, and the granular turbulent phase. In the condensed phase, 
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most of the particles are aggregated around the bottom of the system with an almost close 
packed density and the state is fairly stable; this is because the excitation of the heat bath 
is weak and the potential energy (gravity) is dominant. The collective motion appears on 
the surface and in the bulk of the condensed layer. In the locally fluidized phase, the dense 
layer packing is locally broken by the excitation of the heat bath. The high-speed particles 
are blown upward from the holes, which are the broken layer. However, the location of holes 
is fairly stable. In the granular turbulent phase, the motion of particles becomes very active 
and the kinetic energy is dominant, because the most of particles are highly excited by the 
heat bath. The location of holes becomes unstable in time, and the condensed layer appears 
only temporally. Over the whole region, the flatness f(y = Hf max ) is different from 3, which 
is the value for the Gaussian distribution, but it is remarkable that f(y — Hf max ) becomes 
very large in the locally fluidized phase. The transition from the condensed phase to the 
locally fluidized phase is distinguished by the existence of fluidized holes. On the other hand 
the transition from the locally fluidized phase to the granular turbulent phase is interpreted 
as the destabilization transition of the fluidized hole. 

In the study of the dissipative structure system, such as Rayleigh-Benard convection, it 
is well known that the several phase changes (bifurcations) occur when the external driving 
is increased. However, in dissipative discrete element (granular) systems, there are very few 
studies so far. We believe our study will offer a step for understanding the macroscopic 
universal characters in the studies of fluidization of dissipative (inelastic) discrete element 
system with a heat bath under gravity and constructing the non-equilibrium statistical 
mechanics. 
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FIGURES 

FIG. 1. The maximum packing fraction A max and the excitation ratio fx vs. the driving 
intensity A. The other parameters are fixed at (r, N w , N^) = (0.9, 200, 20). 

FIG. 2. A series of snapshots (interval time is 100r) for the three typical values of A, 
al— a3 (A = 181.8), bl— b3 (A = 250.0), cl— c3 (A = 666.7). The other parameters are fixed 
at (r,N w ,N h ) = (0.9,200,20). 

FIG. 3. The collision rate C to tai/Ttotai y s- the driving intensity A. Exponential fitting curves 
for both LFP and GTP regime are also shown. Each regime has a different exponent. 

FIG. 4. The horizontal velocity PDFs for the three typical As (CP:A = 181.8, LFP:A = 250.0, 
GTP:A = 666.7). The other parameters are fixed at (r,N w ,N h ) = (0.9,200,20). An inset in 
the upper left-hand corner shows the Gaussian nature of the central regions for each distribution 
function. 

FIG. 5. The flatness f(Hf max ) vs. the driving intensity A. The other parameters are fixed 
at (r,N w ,Nh) = (0.9,200,20). The values for Gaussian distribution (/ = 3) and the exponential 
distribution (/ = 6) are indicated in the figure. 

FIG. 6. The velocity PDFs and their fitting functions in LFP and GTP are shown. The 
data are scaled by the standard deviation a. Open circles and filled squares denote that PDFs in 
A = 250.0 with the parameters, xq/cjg = 1.7, (3 = 0.85, and p = 0.25, and A = 476.2 with the 
parameters, xq/gg = 1.7, (3 = 0.85, and p = 0.49, respectively. 

FIG. 7. The values of weight p (the ratio of Gaussian to exponential behavior) for each A are 
shown. 

FIG. 8. The height of maximum values of both packing fraction and flatness for each A are 
shown. 
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FIG. 9. Both maximum packing fractions and packing fractions at the height of Hf max for 
each A in CP are shown. 

FIG. 10. Spatial velocity correlation functions at the height of HAmax in the horizontal 
direction for each A vs. the relative displacement R. 

FIG. 11. The time-dependence of displacement vectors for each particle from the initial 
position are shown at t/r = 100, 200, 300, 400 in A = 25.0. The displacement vectors from the dark 
regions indicate that the particles move toward right-hand side from the initial position. 

FIG. 12. (a) The displacement vectors for each particle from the initial position are shown at 
t/r = 150 in A = 175.4, (b) at t/r = 50 in A = 250.0, and (c) at t/r = 50 in A = 400.0. The 
displacement vectors from the dark regions indicate that the particles move toward right-hand side 
from the initial position. 

FIG. 13. The time-dependence of the horizontal averaged particle current for one particle in 
CP is shown. 

FIG. 14. The spatio-temporal patterns of a surface height of the condensed layer far from the 
critical point ((a) A = 25.0 and (b) 100.0) and near the critical point ((c) A = 188.6) are shown. 
The horizontal and vertical axes correspond to the space expansion (x = ~ L x ) and the time 
evolution (T = ~ 2000r), respectively. 

FIG. 15. The spatio-temporal patterns of the density fluctuations are shown in LFP and GTP. 
(N w = 200; (a) A = 200.0, (b) 250.0, and (c) 454.5). The horizontal and vertical axes correspond 
to the space expansion (x = ~ L x ) and the time evolution (T = ~ 10240r), respectively. The 
black region corresponds to the fluidized holes. 

FIG. 16. The velocity PDFs with fitting function in both N w = 200 and N w = 400 are shown 
(A = 250.0). The data are scaled by the standard deviation a. An inset in the upper right-hand 
corner shows the calculated flatness values as a functions of the weight p for both fitting parameters 



(i.e., (3 = 0.85 and = 0.80) by using Eq. |0. 
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FIG. 17. The spatio-temporal patterns of the density fluctuations are shown in LFP 
(N w = 400; (a) A = 200.0, (b) 250.0, and (c) 333.3 ) and GTP ((d) A = 454.5, (e) 588.2, and (f) 
833.3 ). The horizontal and vertical axes correspond to the space expansion (x = ~ L x ) and the 
time evolution (T = ~ 10240r), respectively. The black region corresponds to the fluidized holes. 
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